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Abstract. - We study the (newtonian) gravitational force distribution arising from a fractal 
set of sources. We show that, in the case of real structures in finite samples, an important role 
is played by morphological properties and finite size effects. For dimensions smaller than d — 1 
(being d the space dimension) the convergence of the net gravitational force is assured by the 
fast decaying of the density, while for fractal dimension D > d — 1 the morphological properties 
of the structure determine the eventual convergence of the force as a function of distance. We 
clarify the role played by the cut-offs of the distribution. Some cosmological implications are 
discussed. 



The aim of the present paper is to discuss the general properties of the gravitational field 
generated by a finite fractal distribution of field sources. This problem is nowadays particularly 
relevant. In fact, there is a general agreement that galaxy distribution exhibits fractal behavior 
up to a certain scale |l], [^]. The eventual presence of a transition scale towards homogeneity 
and the exact value of the fractal dimension are still matters of debate || |5j Qj. Moreover 
it has been observed that cold gas clouds of the interstellar medium have a fractal structure, 
with 1.5<Z3<2ina large range of length scales Q. For this case the general belief is that 
the origin of their fractality lies in turbulence. However recently |tJ it has been pointed out 
that self-gravity itself may be the dominant factor in making clouds fractal. 

Chandrasekhar j^] , has considered the behavior of the newtonian gravitational force proba- 
bility density arising from a poissonian distribution of sources (stars). In this case, in evaluating 
the probability distribution of the force acting on a particle, one supposes that fluctuations are 
subject to the restriction that a constant average density occurs, i.e. the source distribution 
is spatially homogeneous and the density fluctuations obey to the Poisson statistics. Applying 
the Markov's method, it is possible to compute explicitly the force probability density, known 
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as the Holtsmark's distribution. In this case, the probability density of the force F is given by 

W(F) = « (1) 

where F a — (4/15) 2 / 3 (27rGM )n 2 / 3 is the normalizing force (n is the average density of sources, 
M is the mass of each punctual source and G is the gravitational constant), (3 = F/F a is an 
adimensional force and 

2 f°° 

H(j3) = —r / dx exp[-0//?)3] x sin(x) . (2) 

The main result is that, in the thermodynamic limit (V — > oo with rt constant), the force 
distribution has a finite first moment and an infinite variance. This divergence is due to the 
possibility of having two field sources arbitrarily nearby. The approximate solution given by 
the nearest neighbor (n.n.) approximation (i.e. by considering only the effect of the n.n. 
particle) and the exact Holtsmark's distribution (Eq.||) agrees over most of range of F (see 
Fig.l). The region where they differ mostly is when F — > 0. This is due to the fact that a 
weak field arises from a more or less symmetrical distribution of points, and hence the n.n. 
approximation fails. Therefore in the strong field limit we may neglect the contribution to the 
force of far away points, because the main contribution is due to r — > (i.e. it comes from the 
n.n). The root mean square of the force is divergent as in the case of the n.n. approximation, 
and this is due to the fact that the limit r — > is allowed. If there exists a lower cut-off (see 
below) this divergence is erased out. However, due to isotropy, no divergence problem arises 
from the faraway sources even in the limit r — ► oo. 

The derivation of Chandrasekhar cannot be easily extended to the case of fractal distri- 
butions because, in such structures, fluctuations are characterized by long-range correlations. 
A fractal is an instrisic critical structure characterized by self-similarity. The scale-invariant 
properties of such distributions are described by slow power law decreasing density-density 
correlation functions whose exponents are the main characteristic quantities to be studied 
(see H and Q). In this situation an analytic treatment of the force distribution becomes 



very difficult. Vlad 1 1 developed a functional integral approach for evaluating the stochastic 
properties of vectorial additive random fields generated by a variable number of point sources 
obeying to inhomogeneous Poisson statistics. Then he applied these results to the case of 
the gravitational force generated by a fractal distribution of field sources under some strong 
approximations. In particular, in order to compute analytically the force probability density, 
Vlad has not considered the presence of intrinsical fluctuations in the space density. In 
fact, in the computation of the density from a single point one expects to see deviation from 
the average scaling behavior, which are present at any scale (see below). Such a situation 
occurs in any real fractal structure and can be quantified by studying the n-point correlation 
functions [p"i"|| . Instead, the derivation has been done under the assumption that the n + l-point 
correlation function (where the (n + l) th point is the occupied origin) can be written as 
g(x{,X2, "-iXn) oc g{x~i)g{x2)---g{x* n ) , where g(x) is the two-point correlation function. Such 
an approximation is not adequate to describe the effect of morphology and hence in real cases 
the situation is quite different. Instead of EqJ|, Vlad found that the probability density of the 
absolute value F (generalized Holtsmark's distribution) of the field intensity is equal to 



2 f°° 

H(p,D) = — / dx exp[-(x/0)*] 



x sin(x) . (3) 



In this case F D = (4/15) 2/D (2irGM)(DB /(4tt)) 2 / d , where B is a constant characterizing the 
average mass in the unitary sphere and (3 = F/F D . The main change due to the fractal 
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structure is that the scaling exponent in Eq.g is D/2 rather than 3/2. Hence in this case the 
tail of the probability density has a slower decay than in the homogeneous case (see Fig.l). 
This means that the variance of the force is larger for D < 3 than for the D = 3. As we discuss 
below, the case D < 2 is rather well described by Eq.||: this is not the case for D > 2 where 
the n — th points correlations must be taken into account in the case of real structures. An 
important limit is the strong field one (F — > oo). In this case it is possible to show that the 
force distribution of EqJ^ can be reduced to the one derived under the n.n. approximations: 

W nn (F)dF = —(GM)^F-^ exp (-B(GM)% F~%^ dF . (4) 

The n.n. approximation is good for F>Fi, where F\ = (5)"° • GM. From Eq.|], in the 
fractal case, as in the homogeneous one, the divergences of the force moments (the second for 
D < 3 and the first for D < 2) are due only to the fact that two field sources can be arbitrarily 
close. If we impose that a lower cut-off A > (minimal distance between the sources) exists, 
these divergences are erased out. According to the approximations previously mentioned, far 
away points do not produce any divergence neither in the average force nor in its variance. 
This is related to the fact that if the volume involved is large enough with respect to A the 
distribution becomes spherically symmetric in a statistical sense and the anisotropies become 
negligible. 

In the n.n. approximation (Eq.^), we obtain for the average force (F) ~ (1/A) 2 ~ D . Hence 
in the limit A — > 0, (F) is finite for D > 2, while it is infinite for D < 2. The force r.m.s 
value can be computed in the same way and we readily obtain (F 2 ) ~ (1/A) 4 ~ D so that, 
in the limit A — > 0, (F 2 ^ is divergent for any possible value of the fractal dimension such 
that D < 3. Therefore the r.m.s. value of the force strongly depends on the lower cut-off 
A. Consequently, in this context, we expect to see large fluctuations from the average value 
changing the origin point on different masses in the sample. This kind of divergence is due to 
the fact that there is no restriction on the distance Ar between two neighbor field sources, i.e. 
we may have Ar — > 0. To assume A > is quite reasonable for almost every physical problem 
and in particular for the case of galaxy distribution. In fact, if two galaxies were too close, 
they should be subjected to tidal interactions and then should form a binary system. Such a 
situation goes beyond the scope of this paper. 

As we have already mentioned, the Vlad's result has been obtained under an assumption 
which is equivalent to the condition that the mass-length relation (m.l.r.) A^(< r) from a 
single point is exactly a power law function: 

N(< r; j>, fi) = B^r D . (5) 

4"7T 

Eq.^| gives the mass contained in a portion of a sphere of radius r, being a generic fixed point 
at the origin, which covers a solid angle Q in the direction <fi. From this equation we can derive 
directly the number density in the same volume by dividing it for the volume itself (in d = 3) 

n(<r;ln) = ^-r D - 3 (6) 

It is important to notice, instead, that the m.l.r. has a genuine power law behavior only when it 
is averaged over all the points x of the structure, obtaining the so-called two-point correlation 
function (tpcf). From a single point there can be important intrinsical fluctuations around the 
power law behavior JTT[ . In fact, since a fractal is an intrinsically critical system with strong 
correlations on all scales, we expect that the density from one point shows large fluctuations, 
with respect to the tpcf, changing r with fixed q> or vice-versa. In general these intrinsic 
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fluctuations can be characterized as log-periodic corrections to scaling [|Tl|] . These fluctuations 
can be present at any scale (i.e. whatever is r) and they can be angularly correlated, that is 
the quantity 

dn(\$i - h\;r) = ((n(< r-J u 6Q) - n(< r;0 2 ,^)) 2 ) (7) 

depends on the angular distance between <f>i and <j>2- As we show below, these fluctuations 
can play a crucial role in the determination of the force distribution. In order to consider, 
in a simple way, the important effect of morphology, we study the r.m.s force for a generic 
random fractal distribution in the three dimensional Euclidean Space. We consider explicitly 
the existence of a minimal distance (i.e. lower cut-off) A and a maximal distance (i.e. upper 
cut-off) R s from the origin, and then we derive the asymptotic behavior for A — > and 
R s — > oo. Let F be the total force acting on a particle due to all the other field sources having 
the same unitary mass (GM — I). We consider the contribution to the force by all the field 
sources in the spherical shell [A, R s ] 



sr 



F{^Rs)= ^ $ (8) 

A<\fi\<R s 1 

where N is the total number of the sources in the shell considered and fi is the position of 
the i — th source with respect to the one in the origin. Under the assumption that the mass 
distribution is statistically isotropic, we have 

lim (f(A,R s )) =0. (9) 
The statistical average of the force square modulus is given by: 

'\F(A,R s )\*) = (F(A,R s f) = ( £ £) + ( E CJ W) ( 10 ) 

where 8ij is the angle formed by fi and fj . Here the average < ... > is made over an ensemble 
of different realizations of the sample or, when this is not possible, over all the possible origins 
in a single sample. Obviously, the first term is simple to analyze, while the second one deserves 
a more careful discussion because it involves the three-point correlation function. The first 
term can be written as 

/ ^ l\ f R ' , r D ~ l DB I / (A\ 4 - D \ , , 

In Eqjl| we have considered that < ... > leaves only the power law behavior in the m.l.r. In 
fact we can write X^i=i pt = Ia r d 3 rn(f)/r 4 , where n(f) = ^2^ =1 S(f — fi) is the number 
density, which can be different from a power law and show persistent fluctuations around its 
average behavior, and < ... > erases these fluctuations as well as the dependence from the 
direction of f. Analogously, the second term in Eq.[[o| can be written as 

/ l ^^' N cos(%) \ f [ ^,^cos(6ii2) 

\A<|ri|<_R s * 3 I 



d'nd^pin,^)—^ (12) 

VJV r l r 2 



where p(fi, r*2)d 3 rid 3 r2 is equal to N(N — 1) times the probability to have a particle in 
the infinitesimal volume d 3 r± and, at the same time, another in the volume d?T2, i.e. the 
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number of pairs in which the first particle is in d?r\ and the second in d 3 r 2 . The function 
vi^ii^i) is strictly related to the three-point correlation function in which the first point 
is the origin and the others are r*i and r 2 - Such a function is usually very difficult to 
be computed for a generic fractal structure [ ftT| . Note that p(fi,r 2 ) ~ (n(fi,r 2 )) where 
n(fi, r 2 ) = J2i j ^(^ — r i)^( r j — • We expect that p(ri,f 2 ) depends only on n and r 2 and 
the angle 9 between the two directions because of the average < ... >. We make the following 
approximation, by using the tpcf and the angular correlation function 

/ p(fi > r 2 )dClxdn 2 ~ (4tt) 2 (n(n)> (n(r 2 )) V [0, — ) d9 (13) 

^6»<|£li — £l 2 |<f+de V r 2/ 

where T(6,r\/r 2 ) is the probability that the angle between fi and f 2 is in the range [9,8 + d8], 
conditioned to the fact that the distances of the two sources are r\ and r 2 . The dependence 
of T{8,r\/r 2 ) on the ratio r±/r 2 is in general a peculiar property of the fractal studied. This 
problem has been studied in detail by f0] |. In what follows, we neglect the dependence on 
the ratio r\jr 2 , which can be treated as a perturbation ]TT[ . This approximation consists in 
substituting T(8,rx/r 2 ) with the angular correlation function T(8) normalized to unity. T(8) 
is obtained by the analysis of the angular projection of the sample. 
Under this approximation we can write Eq.|l^ as 

/ cos(%) \ {DB fC(D) 1 ( ( A J A ^ 2 - DN 



\ ^ r 2 r 2 l~ (2-D) 2 A*-"> \ 1+ \R a ) 2 UJ J (W) 

\A<\fi\<R s % J / V ' \ \ »/ \ »/ / 

where C(D) = £ T(9) cos (9)d9 . The presence of fluctuations in the m.l.r. from one point 
are crucial in determining the behavior of T(6) and then the convergence properties of (f 2 )- 
To clarify this point, consider the quantity dn{\4>\ — 4> 2 \]t) in EqJ^. It is simple to show 
analytically that this quantity is strictly related to T{9) where 6 = \(f>i — (f> 2 \ is the angular 
distance between </>i and <f>2- At this point is convenient to distinguish the cases D < 2 and 
D > 2. For D < 2, even in absence of fluctuations in the m.l.r., we have C(D) > (the 
effect of fluctuations is only a change in the non-vanishing value of C(D)), while for D > 2, 
C(D) > only because of these fluctuations. In the latter case and in absence of fluctuations, 
one would have a case very similar to the homogeneous one |^] in which T(9) is constant and 
therefore C(D) =0. 

Consequently, for D < 2, the leading behaviors of (F 2 } are: 

(F 2 ) ~ a{B,D)A-( 4 - D ^ /orA->0 

and R s finite 

(F 2 )cFi-2l3(B,D)(±) 2 ~ D /orf - oo (15) 

and A > 

In Eq.|l| a(B,D) = BD/(4 - D) , 0(B,D) = {DB) 2 C{D)/[{2 - D) 2 A 4 - 2D and F^ = 
a(B 1 D)A-^- D '> + (3{B,D) is the finite asymptotic value of (F 2 ) for R s /A oo. Obviously 
for intermediate ranges of R s /A, other powers of it in Eqjll] and EqJIJ are important. 

The case D > 2 is very different and the effect of eventual fluctuations in the m.l.r. from 
one point can be crucial for the convergence of (F 2 ) in the limit R s / A — > oo. As a first step 
let us suppose that the m.l.r. from one point is a pure power law. In this case the function 
T(9) is expected to be a constant and then C{D) — 0. As a consequence we have that the 
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second term in Eq.^ is equal to zero and (-F 2 ) is simply given by 

This quantity diverges for A — > as A - ' 4- - ' and for A > it converges for R s /A — > oo. 

In the case in which fluctuations around the power-law relation in the m.l.r. are present 
(and this is the more general case), the behavior is more complex. In this case T(6) is not 
anymore constant, since we expect that the effect of these fluctuations is to introduce some 
angular correlations at least at small angles. The reason is that these fluctuations are due to 
the sequence of voids and clusters of masses in a fractal sample. These sequences are not the 
same in different directions, but they are angularly correlated and the average < ... > does 
not erase these correlations. In such a case (F 2 ) is given by 



(F 2 ) = 



DB I ( ( A\ 4 ~ D \ (DB) 2 C(D) 1 ( ( R s N 



4-DA 4 - D \ \R S J r (D-2) 2 A 4 - 2D \ \A 




From Eq.(p7[) we can see that, for R s finite, (-F 2 ) diverges as A~( 4 ~-°) (as in the case D < 2) 
when A — ► 0, and for A > it diverges as (R s / A) 2D ~ 4 in the limit R s /A — + oo. This last 
behavior is in striking contrast with the Vlad's one for the same values of D; this is due to 
the fact that as mentioned Vlad does not consider the fluctuations in the m.l.r. from a single 
point. At intermediate scales the behavior is more complex and depend on the values of the 
parameter D, B and the value of C(D). 

We have tested the results obtained in Eq.20 by numerical simulations. We have generated 
artificial fractals by means of different algorithms: random cantor set, levy flight and random 
walk M. With these algorithms the fractal dimension can be tuned appropriately (i.e. the 
tpcf), while the morphological properties of the structure (and hence the higher order correla- 
tions) are determined by the kind of iterative process choosen. For this reason the simulations 
do not have the same morphology of real galaxy distribution. A more realistic study for 
the case of galaxy structures is in progress, and here we are concerned with a more general 
problem. The agreement with Eq.20 is rather good for both cases D < 2 and D > 2. We stress 
again that the divergence of Eq.20 as a function of R s is strictly related to the morphological 
properties of the structure considered. 

The clarification of the properties of the gravitational field in fractal structures is partic- 
ularly relevant for the studies of the matter distribution in the universe. It is commonly 
believed that the peculiar velocity field (derived as a distortion from a pure linear Hubble 
law) is generated by local matter inhomogeneities in the distribution of galaxies |0|. A 
standard hypothesis in this framework is that, on large enough scale, the linear perturbation 
theory holds, as density fluctuations are small enought, and by using such an approximation 
one may "reconstruct" the peculiar velocity field from the matter one or vice-versa. Clearly, 
if matter distribution shows fractal properties, the linear theory fails. A related problem 
consists in the computation of the net gravitational force, due to the local inhomogeneities 
of matter distribution, acting on our galaxy (the dipole) [jl2| . Here we have shown that the 
eventual convergence or divergence of the dipole is related to subtle morphological properties 
of the sources distribution and not only to the prtoperties of the tpcf. In the universe we 
observe an irregular distribution of galaxies characterized by fractal properties , and hence 
its relationship with the peculiar velocity field must be much more complex than the one 
expected in the linear theory p^[ . Such a situation deserves further studies, and in this paper 
we have considered some basic problems related to it. 
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Fig. 1. - (Upper panel) Force distribution due the n.n. (crosses) and due to all the field sources 
(diamonds) for an homogeneous sample with n = N/V = 2.39 ■ 10 4 in d = 3. The dotted line 
represents the force distribution W(F) computed by the n.n. approximation, while the dashed line 
is the Holtsmark's distribution (Eq.2). The agreement is quite satisfactory at strong field, and there 
is a deviation at weak field. (Bottom Panel) The same for a fractal with D = 1.98 generated by a 
random walk in d = 3. In this case we have used the generalized Holzmark distribution and its n.n. 
approximation for the fits (dotted and dashed lines) 
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Fig. 2. - (Left Panel) Behavior of ^/(F(A,R a ) 2 ) as a function of R a for a fractal with D = 1.5 
generated by the random cantor set algorithm. The dashed line represents the first term in Eq.20, 
while the solid line both contributions in Eq.20. The agreement is rather good. (Right Panel) The 
same for a fractal with D = 2.7. The effect of the angular correlations is now important and determine 
the divergence of y {F(A, Rs) 2 }. The error bars have been computed over an average of 5 different 
realization of the same structure. 



